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1 .  Introduction 

The  objective  of  this  project  was  to  investigate  the  feasibility  of  imaging  subsurface 
acoustic  sources  generated  artificially  within  an  underground  facility.  One  application 
identified  as  important  by  the  Air  Force  is  source  location  within  a  buried  facility  based 
on  the  passive  monitoring  of  acoustic  emissions  from  both  stationary  and  moving 
equipment  within  such  a  facility.  Another  potential  objective  is  to  identify  the  type  of 
source  on  the  basis  of  its  noise  spectrum,  although  source  characterization  was  beyond 
the  scope  of  this  project;  rather,  our  goal  was  to  examine  and  demonstrate  the  possibility 
of  source  imaging.  The  imaging  techniques  studied  are,  in  principle,  capable  of 
reconstructing  an  intensity  map,  or  image,  of  an  extended  source  distribution  using 
signals  passively  recorded  by  a  geophone  array  deployed  above  ground. 

Several  approaches  to  passive  seismic  imaging  have  been  studied  in  recent  years.  A 
method  investigated  by  Shuster  et  al.  [1]  and  Wapenaar  et  al.  [2-5]  termed  "seismic 
interferometry,"  exploits  the  interference  of  a  direct  arrival  with  a  scattered  or  secondary 
arrival,  which  may,  for  example,  originate  from  baekscattering  from  the  surface.  This 
approach  is  an  elaboration  of  an  idea  by  Claerbout  [6]  relating  transmission  and  reflection 
seismograms  using  autocorrelated  data  in  which  one  component  originates  from  a  free 
surface  reflection.  If  the  subsurface  source  is  impulsive  or  very  broadband,  then  the 
cross-correlation  between  the  direct  and  secondary  arrivals  can  bear  a  strong  resemblance 
to  a  reflection  seismogram. 

An  essential  ingredient  in  the  above  methods  is  the  interference  between  a  direct  and  one 
or  more  secondaiy  arrivals.  The  passive  imaging  methods  investigated  in  the  current 
project,  however,  differ  in  that  they  employ  only  the  direct  arrivals;  no  interference  with  a 
"reference  wave"  is  required.  In  particular,  two  imaging  algorithms  of  this  kind  have  been 
previously  published  [7,8]  by  the  PI  and  one  outcome  of  this  work  was  a  clearer 
understanding  of  the  relationship  between  these  two  techniques  [9].  The  earliest  of  these 
algorithms  was  derived  for  imaging  acoustic  emission  sources  in  materials  [7],  while  the 
second  algorithm  was  intended  for  seismic  imaging  [8].  It  is  important  to  note  that  these 
algorithms  do  not  use  time  delays  across  an  array  to  "triangulate"  on  the  source,  which  is 
inherently  a  nonlinear  procedure.  Rather,  the  algorithms  are  truly  imaging  methods, 
capable  of  mapping  a  spatially  extended  source.  The  imaging  algorithms  are  linear  in  the 
intensity,  and  thus  the  source  need  not  be  spatially  localized.  In  order  for  the  algorithms 
to  be  linear,  the  assumption  of  a  spatially  incoherent  source  is  required,  meaning  that  the 
signals  emitted  from  spatially  distinct  points  within  an  extended  source  are  assumed  to  be 
statistically  uncorrelated. 


In  one  imaging  algorithm  [7],  termed  "time-exposure  acoustics,"  passive  data  acquired  by 
an  array  of  geophones  are  used  to  reconstruct  a  source  image  by  backpropagating  the 
signals  into  image  space.  This  procedure  bears  some  resemblance  to  time-domain 
migration.  But,  in  contrast  to  the  active  source  problem  in  which  signals  are  transmitted 
into  the  ground,  in  the  passive  problem  no  time  origin  exists  with  respect  to  which  the 
propagation  delays  can  be  referenced.  That  is,  the  signals  are  assumed  emitted  from  the 
source  in  a  continuous  and  random  fashion.  An  essential  idea  described  in  [7]  is  to 
average  with  respect  to  this  time  origin.  In  particular,  multiple  images  are  reconstructed 
each  with  a  different  assumed  time  origin  (the  time  at  which  the  signal  is  emitted)  and 
then  these  images  are  averaged  on  an  intensity  basis.  This  approach  exploits  the  fact  that 
the  signal,  either  emitted  from  or  scattered  by  a  point  in  the  subsurface,  will  create  a 
wavefront  that  will  manifest  itself  as  a  correlated  signal  across  a  receiving  aperture  on  the 
surface.  By  coherently  summing  this  signal  and  backpropagating  to  multiple  source 
locations,  one  can  develop  an  image  of  source  intensity.  As  discussed  in  [7]  image 
resolution  is  determined  essentially  by  two  factors:  the  source  bandwidth  and  the  array 
size. 

The  above  method  has  been  referred  to  as  time-exposure  acoustics  (TEA)  since  the 
source  is  continuous  and  ever  present,  and,  as  a  result,  an  image  can  be  built  up  over  time 
as  data  are  continually  accumulated.  The  underlying  idea  is  to  exploit  the  spatial 
coherence  of  signals  which  are  correlated  across  an  array  to  construct  an  image.  When 
properly  averaged,  random  fluctuations  in  an  ensemble  of  images  will  diminish  over  time 
as  new  data  are  incorporated. 

In  the  second  source  reconstruction  algorithm  [8]  the  source  imaging  problem  is 
formulated  as  a  tomographic  reconstruction  problem.  Here,  signals  recorded  by  different 
pairs  of  geophones  are  cross  correlated.  The  resultant  correlation  peak  consists  of  an 
integration  of  the  source  intensity  lying  along  a  hyperbolic  path  (in  a  2-D  problem)  or  on 
the  surface  of  a  hyperboloid  (in  a  3-D  problem)  in  which  the  coordinates  of  the  two 
geophones  define  the  foci  of  the  hyperbola  or  hyperboloid.  In  this  formulation,  the 
correlation  peaks  may  be  regarded  as  the  line  integral  of  source  intensity  over  this 
hyperbolic  path  (in  2-D)  or  as  the  surface  integral  of  source  intensity  over  a  hyperboloid 
(in  3-D).  It  is  shown  in  [9]  that  this  method  can  be  manipulated  into  a  form  that  closely 
resembles  the  TEA  algorithm,  although  this  relationship  may  not  be  immediately 
obvious. 

In  the  derivation  of  both  algorithms,  we  assumed  a  scalar  wave  (acoustic)  problem  and  a 
constant  velocity;  also  mode  conversion  effects  at  the  Earth’s  surface  were  neglected. 
These  assumptions  are  made  for  simplicity,  but  the  essential  ideas  presented  here  can  be 
generalized  in  principle  to  more  complex  (and  realistic)  earth  models. 

In  this  final  report,  a  review  of  the  fundamental  imaging  equations  of  both  algorithms  is 
given  in  the  Appendix.  A  detailed  derivation  of  these  equations  and  how  the  two 
algorithms  are  related  can  be  found  in  [9].  The  latter  article  was  one  outcome  of  this 
research.  We  shall  also  review  below  some  key  images  derived  from  field  data.  We 
finally  present  new  results  obtained  from  data  recorded  over  a  coalmine  during  mining 
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operations.  These  last  results  demonstrate  that  actual  images  can  be  derived  under 
realistic  conditions.  We  conclude  with  suggestions  for  new  areas  of  investigation  that 
may  further  improve  image  quality. 

2.  Images  of  Points  Sources  —  Experimental  Setup 

Below,  we  compare  the  two  reconstruction  algorithms  by  applying  them  to  the  same  field 
data.  The  data  in  these  tests  were  acquired  with  a  geophone  array  that  recorded  surface 
waves  generated  by  a  hammer  source  at  a  single  location.  Since  the  source  and  array  are 
on  the  surface,  the  images  are  derived  from  surface  (Rayleigh)  waves  only.  These  waves 
are  dispersive  and  thus  not  in  principle  optimal  for  imaging  since  the  algorithms  are  not 
designed  to  account  for  wave  dispersion.  Subsurface  body  waves  in  3-D  are  not 
dispersive  and  should,  in  theory,  provide  better  images.  Nevertheless,  as  shown  below, 
quite  good  images  of  the  point  sources  were  obtained  using  surface  waves. 

As  a  first  test,  seismic  data  were  recorded  by  a  seismograph  attached  to  an  8  geophone 
array  with  a  15  foot  spacing  between  elements.  The  seismograph  was  not  triggered,  but 
merely  turned  on  for  30  seconds  at  a  time;  a  large  number  of  these  30  second  intervals 
were  recorded.  The  source  was  a  hammer  blow  on  a  metal  plate,  which  created  surface 
waves  that  impinged  upon  the  array.  The  source  location  was  45  feet  from  the  array.  A 
program  was  written  to  convert  the  SEG-2  file  format  generated  by  the  seismograph  into 
a  binary  file  capable  of  being  read  by  a  MatLab  program,  which  performed  the  image 
reconstructions. 

Figure  1  shows  the  raw  seismic  data  generated  by  the  hammer  source  recorded  by  the  8 
element  array.  The  power  spectrum  of  the  signal  recorded  by  one  geophone  is  shown  in 
Figure  2,  which  indicates  a  receiver  bandwidth  of  about  5  Hz  to  50  Hz.  The  surface  sound 
velocity  was  estimated  from  the  curvature  of  the  impinging  wavefront  and  knowledge  of 
the  source  location  to  be  approximately  600  ft/sec.  In  practice,  when  the  source  location 
is  unknown,  the  velocity  cannot  be  measured  directly.  A  simple  approach  to  estimating 
the  correct  velocity  is  to  reconstruct  multiple  images  using  a  range  of  velocities  and  then 
to  select  the  velocity  that  provides  the  sharpest  image.  Various  measures  of  image 
sharpness  can  be  used  for  this  purpose,  e.g.,  borrowed  from  the  field  of  adaptive  optics. 
Sharpness  criteria,  for  example,  can  be  based  on  the  mean  entropy  of  the  image  or  the 
mean  spatial-frequency  bandwidth  of  the  image. 
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Figure  J .  Data  acquired  over  a  30  second  time  window  with  8  geophones.  Over  this 
period  of  time,  the  ground  was  struck  with  a  hammer  12  times.  No  seismograph 
triggering  was  used  in  acquiring  the  data.  The  sampling  interval  is  4  msec. 


Spectrum  of  Received  Signal 


Figure  2.  The  power  spectrum  of  a  30  second  time  record  recorded  by  one  geophone. 
This  shows  that  the  bandwidth  is  roughly  5  Hz  to  50  Hz. 
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3.  Reconstructed  Images  from  Field  Data 

Both  the  time-domain  migration  and  the  cross-correlation  algorithms  were  applied  to  two 
data  sets  (set  1  and  set  2)  recorded  by  the  array.  Figures  3  is  a  reconstructed  image  of  the 
hammer  source  using  the  time-domain  migration  algorithm  using  data  set  1.  Figure  4  is 
the  image  of  the  source  using  the  same  data  set,  but  employing  the  cross-correlation 
algorithm.  In  these  reconstructions,  30  seconds  of  data  were  used.  The  cross-correlation 
images  are  similar  to  the  time-domain  migration  images  with  the  latter  (Fig.  3)  appearing 
to  have  somewhat  better  definition. 

Figures  5  and  6  are  images  using  data  set  2  with  the  same  two  algorithms,  respectively. 
This  data  set  had  somewhat  superior  signal  to  noise,  and  as  a  result  the  quality  of  the 
images  appears  somewfrat  better  than  those  derived  from  data  set  1 . 

Note  that,  in  all  of  these  images,  the  longitudinal  resolution  (resolution  in  the  direction 
perpendicular  to  the  army)  is  elongated  compared  to  the  lateral  resolution.  This  is 
expected,  since  longitudinal  resolution  is  determined  solely  by  our  ability  to  measure  the 
curvature  of  the  impinging  wavefront.  Lateral  resolution,  on  the  other  hand,  is  primarily 
determined  by  the  bandwidth  of  the  noise  and  the  array  length.  Resolution  can  be 
improved  if  one  uses  two  linear  arrays  perpendicular  to  each  other.  This  was  done  in  a 
second  experiment  in  which  two  arrays,  each  with  seven  elements,  were  arranged  in  the 
shape  of  an  ML."  The  geophones  were  again  spaced  at  15  foot  intervals.  The  image  of  a 
hammer  source  45  ft.  away  from  both  arrays  is  shown  in  Figure  7  in  which  the  time- 
domain  migration  algorithm  was  used.  Using  the  same  data,  the  image  reconstructed 
using  the  cross-correlation  algorithm  is  shown  in  Figure  8.  With  the  L-shaped  array,  note 
that  the  resolution  is  tighter  and  more  symmetrical  compared  to  the  images  of  Figure  3 
through  6.  Figure  9  show  images  using  the  time-domain  migration  algorithm  with  data 
from  only  a  single  leg  of  the  array.  In  Figure  9,  the  array  is  in  the  y-direction  and  the 
image  is  elongated  in  the  direction  perpendicular  to  the  array  as  expected.  Figure  10 
shows  images  using  the  time-domain  migration  algorithm  with  data  from  the  other  leg  of 
the  array  (oriented  in  the  x-direction).  Again,  elongation  occurs  perpendicular  to  the 
array.  When  the  data  from  both  legs  are  combined,  we  obtain  the  symmetrical  image 
shown  in  Figure  7. 


5 


Figure  3.  Image  obtained  using  the  time-domain  migration  algorithm  (data  set  1). 


Figure  4.  Contour  image  obtained  using  the  cross-correlation  algorithm  (data  set  1). 
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Figure  5.  Image  obtained  using  the  time-domain  migration  algorithm  (data  set  2). 


Figure  6.  Image  obtained  using  the  cross-correlation  algorithm  (data  set  2). 
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Figure  7.  Image  obtained  from  the  L-shaped  array  data  using  the  time-domain  migration 
algorithm. 


Figure  8.  Image  obtained  from  the  L-shaped  array  data  using  the  cross-correlation 
algorithm. 
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Figure  9.  Image  obtained  from  the  first  leg  (7  geophones  in  the  y  direction)  of  the 
L-shaped  array  using  the  time-domain  migration  algorithm. 


0  20  40  60  80 


Figure  10.  Image  obtained  from  the  second  leg  (7  geophones  in  the  x  direction)  of  the 
L-shaped  array  using  the  time-domain  migration  algorithm. 
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4.  Coalmine  Images 

Two  orientations  of  a  linear  array  of  geophones  were  used  to  record  seismic  data 
generated  by  mining  operations  over  a  coalmine  in  Salineville,  Ohio  on  July  28,  2005. 

The  arrays  were  laid  out  perpendicular  to  each  other  and  an  attempt  was  made  to  center 
the  arrays  directly  over  the  region  of  the  mine  undergoing  the  cutting  operation.  The 
active  mining  face  (the  source  region)  was  350  ft  below  the  surface.  The  array  length  was 
420  ft  with  29  geophones  spaced  at  1 5  ft  intervals.  Data  were  recorded  with  the  first  array 
(Line  1)  over  a  period  of  about  40  minutes.  For  Lined,  40  30-second  data  records  were 
recorded.  The  second  array  orientation  (Line  2)  was  perpendicular  to  Line  1  with  the 
same  length  and  number  of  geophones.  Line  2  recorded  38  30-second  data  records.  Line 
1  was  oriented  parallel  to  the  mine  tunnel,  while  Line  2  was  oriented  perpendicular  to  the 
mine  tunnel.  A  sampling  interval  of  2  ms  was  used.  Figure  1 1  shows  the  traces  recorded 
by  the  29  geophones  from  one  30  second  data  record.  Notice  the  rather  sudden  onset  of 
activity  toward  the  end  of  this  record.  Figure  12  is  the  spectrum  of  a  single  trace, 
indicating  some  energy  extending  beyond  100  Hz,  although  most  of  the  signal  energy  is 
confined  to  below  50  Hz. 

Both  imaging  algorithms  yielded  very  similar  results,  so  we  show  below  images  obtained 
using  the  time-domain  migration  algorithm  only.  Figures  13a  and  13b  show  40  images 
reconstructed  using  the  Line  1  data.  Figure  13c  shows  the  image  obtained  by  combining 
all  40  records  of  data.  Similarly,  Figures  14a  and  1 4b  show  38  images  reconstructed 
using  Line  2.  Figure  14c  shows  the  image  obtained  by  combining  all  38  records  of  the 
Line  2  data.  Some  variability  in  the  images  shown  in  Figures  13a,  13b,  14a  and  14b  was 
to  be  expected  since  the  operation  of  the  mining  machinery  was  intermittent  over  the  40 
minute  recording  interval. 

Using  the  Line  1  data  set,  it  was  found  that  a  sound  velocity  of  9000  ft/sec  placed  the 
source  image  at  the  correct  depth  (350  ft).  This  is  within  the  range  of  measured  sound 
speeds  for  coal,  estimated  to  be  from  about  6000  to  1  1 ,000  ft/sec  [10].  Flowever,  for  the 
Line  2  data  set,  a  higher  velocity  of  12,000  ft/sec  was  required  to  place  the  source  at  the 
correct  depth.  This  difference  in  apparent  velocity  may  be  due  to  several  factors.  One 
possibility  is  the  presence  of  higher  velocity  material  (e.g.,  sandstone)  overlying  the  coal 
deposits.  Another  may  be  due  to  the  fact  that  the  Line  1  data  were  recorded  with  the  array 
parallel  to  the  mine  tunnel,  which  could  have  artificially  reduced  the  apparent  velocity 
due  to  interaction  with  the  tunnel  (e.g.,  the  tunnel  acting  as  a  wave  guide  creating  multi- 
path  effects).  Future  calibration  tests  using  a  known  source  (e.g.,  a  hammer  blow)  at  a 
known  location  within  the  mine  could  help  explain  the  reasons  for  this  velocity 
difference. 
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Figure  1 J .  Traces  of  the  raw  data  recorded  by  29  geophones  over  30  seconds.  Notice  the 
onset  of  activity  toward  the  end  of  this  record. 
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Figure  12.  Spectrum  of  a  single  trace. 
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Figure  13a.  Each  image  shown  was  derived  from  a  separate  30  second  data  record  from 
the  Line  1  data  set.  Depth  is  downward.  The  peaks  of  the  images  are  consistent  with  a 
localized  source  at  a  depth  of  350  ft.  The  variability  in  the  images  is  a  result  of  the 
mining  activity  being  intermittent. 
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Figure  13b.  20  more  images  derived  from  the  Line  1  data  set. 
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combined  data 
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Figure  ]  3c.  The  image  obtained  by  combining  all  40  data  sets  from  the  Line  1  array 
orientation.  The  source  was  at  a  depth  of  350  feet. 
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Figure  14a.  Each  image  shown  was  derived  from  a  separate  30  second  data  record  from 
the  Line  2  data  set. 
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Figure  14b.  18  more  images  derived  from  the  Line  2  data  set. 
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Figure  14c.  The  image  obtained  by  combining  all  38  data  sets  from  the  Line  2  array 
orientation. 
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5.  Recommendations  and  Conclusions 


The  coalmine  study  demonstrates  that  passive  imaging  is  feasible  under  realistic 
conditions,  but  improvements  are  certainly  possible.  For  example,  some  noticeable  phase 
rotation  of  the  seismic  traces  across  the  array  is  evident  when  the  traces  are  plotted  (e.g., 
as  in  Figure  1 1  but  at  higher  magnification).  The  phase  rotation  is  probably  due  in  part  to 
dispersion  of  the  seismic  waves  as  a  function  of  path  length  (greater  path  lengths  imply 
greater  dispersion).  Phase  rotation  will  degrade  the  resolution  of  the  image.  Simulations 
show  that  images  derived  from  idealized  data  in  the  absence  of  phase  rotation  are  of 
higher  quality.  No  phase  rotation  implies  that  the  signals  recorded  at  different  geophones 
are  merely  delayed  versions  of  one  another,  where  the  delay  is  proportional  to  path  length 
differences  from  the  source  to  the  receivers.  Phase  rotation  due  to  dispersion  can  be 
corrected  to  some  extend  by  inverse  filtering,  but  was  not  attempted  here.  The  reason  for 
\  this  is  that  rotation  due  to  pure  dispersion  due  only  to  path  length  differences  was  not 
clearly  evident  in  the  raw  data,  since  systematic  phase  distortion  as  a  function  of  path 
length  was  not  seen  across  the  array.  This  seems  to  imply  that  the  phase  distortion  is 
probably  due  primarily  to  velocity’  inhomogeneities  along  the  different  propagation  paths. 
Partial  compensation  for  this  might  be  achieved  by  incorporating  a  layered-earth  model 
|  into  the  imaging  algorithms,  although  this  would  require  some  a  priori  knowledge  of  the 
local  stratigraphy. 

Another  improvement  that  could  be  made  in  future  measurements  would  be  to  reduce  the 
sampling  interval  to  well  below  the  2  ms  used  here,  since  the  difference  in  arrival  time 
between  adjacent  phones  was  typically  smaller  than  2  ms  (but  sufficient  to  resolve  delay 
differences  between  non-adjaeent  phones). 

Since  the  arrays  were  linear,  only  2-D  images  were  obtained  with  no  spatial  resolution 
perpendicular  to  the  array.  Fully  3-D  resolution  would  require  a  2-D  surface  array.  The 
optimum  way  to  distribute  a  limited  number  of  geophones  in  2-D  to  achieve  the  best 
resolution  would  also  require  further  study. 

Another  area  that  could  be  investigated  to  improve  image  quality  and  resolution  is  the  use 
of  adaptive  imaging  techniques,  similar  in  concept  to  adaptive  techniques  used  in 
astronomy  to  compensate  for  random  phase  eirors  arising  from  propagation  through  an 
inhomogeneous  atmosphere  [11-12].  This  approach  works  well  when  one  can  identify  a 
single  isolated  source  whose  image  can  be  sharpened  through  an  adaptive  process 
involving  the  application  of  phase  delays  to  the  signals  recorded  over  an  aperture  (in  our 
case,  the  sensor  array).  These  phase  delays  are  selected  to  maximize  some  measure  of 
image  sharpness.  Any  sources  that  are  nearby  the  "reference"  source  will  also  undergo 
sharpening  as  long  as  the  propagation  paths  from  nearby  sources  are  similarly  affected  by 
local  inhomogeneities  in  the  propagating  medium.  This  region  for  which  sharpening 
occurs  is  known  as  an  isoplanatic  patch  in  optics.  As  far  as  we  are  aware,  adaptive 
techniques  have  not  been  applied  to  passive  or  active  seismic  imaging.  Certainly,  the 
conditions  are  not  as  favorable  as  in  astronomy  (in  which  a  single  bright  star  typically 
serves  as  the  isolated  reference  source),  but  if  a  dominant  and  isolated  seismic  source  can 
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Appendix:  Passive  Imaging  Algorithms 

A.  Time-Domain  Migration  Algorithm 

Letting  u{r„  I !  denote  the  signal  recorded  by  the  n-th  geophone,  the  imaging  algorithm 
is  given  by  the  formula 


"I 


where 


«„(r)  =  |r  —  rr,|  u(r„,  |r  —  r„|/c  l  /0). 


is  the  backprojected  signal  and  the  angle  brackets  denote  a  time  average  with  respect  to 


B .  Cross-Correlation  Algorithm 

Define  the  cross-correlation  between  the  signals  recorded  by  the  m-th  and  n-th  geophone 
as 


rm,(r„.  r t)u(r„ „  /'))fo 


Then  the  cross-correlation  image  reconstruction  algorithm  is  given  by 
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